/*!  \file  SteVX.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      07/07/2013 12:37:13 PM
 *   \copyright GNU Public License.
 */

#include "SteVX.h"
#include "Matrix.h"
#include "setThreads.h"
extern "C" {
void dstevx_( 
const char   *  jobz, 
const char   * range, 
const int    *     n, 
      double *     d, 
      double *     e, 
const double * vl, 
const double * vu, 
const int    * il, 
const int    * iu, 
const double * abstol,
int    * m, 
double * w, 
double * z, 
const int    * ldz, 
double * work, 
int    * iwork, 
int    * ifail,
int    * info );
}
SteVX::SteVX(): TriEigen("Lapack dstevx") {	
	setCompileTime(__DATE__, __TIME__); 
	canBeParallel = false;
}
bool SteVX::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
#include <iostream>
using namespace std;
void SteVX::compute(const Matrix & T, Matrix & Sigma, Matrix & Q ){
	reset();
	assert(T.isTriDiagonal() && Sigma.isDiagonal() && Q.isOrthogonal());
	char jobz = 'V', range = 'A'; int n = T.cols(), info = -1, ldq = Q.lda();
	double * e = NULL; allocate(e, n); 
	double * d = NULL; allocate(d, n); 
	memcpy(d, T.getData()    ,    n * sizeof(*d));
	memcpy(e, T.getOffDiagonal(), (n-1)* sizeof(*d));
	double vl = 0; double vu = 1; 
	int il = 1; int iu = n; int m;
	double abstol = 0;  
	double * w = Sigma.getModData();
	double * q = Q.getModData(); Q.makeColMajor();
	int lwork = 5 * n, liwork = 5 * n;
	success = false;
	double *  work = NULL; allocate( work, lwork);
	int    * iwork = NULL; allocate(iwork, liwork);
	int    * ifail = NULL; allocate(ifail, n);
	dstevx_(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &abstol, &m, w, q, 
			&ldq, work, iwork, ifail,	&info);
	success = info == 0;
	delocate(ifail);
	delocate( work);
	delocate(iwork);
	delocate(d); 
	delocate(e); 
};
BiNESteVX::BiNESteVX() : BiSVD( "NE dstevx" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteVX()));
	if (eigen) setCompileTime(eigen->getCompileTime());
};
BiNE2SteVX::BiNE2SteVX() : BiSVD( "NE2 dstevx" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteVX()));
	setCompileTime(__DATE__, __TIME__);
};
//REGISTER(SteVX);
//REGISTER(BiNESteVX);
//REGISTER(BiNE2SteVX);
